Theory of the random potential and conductivity at the surface of a topological 

insulator 



Brian Skinner and B. I. Shklovskii 

Fine Theoretical Physics Institute, University of Minnesota, Minneapolis, MN 55455, USA 

(Dated: January 22, 2013) 

We study the disorder potential induced by random Coulomb impurities at the surface of a 
topological insulator (TI). We use a simple model in which positive and negative impurities are 
distributed uniformly throughout the bulk of the TI, and we derive the magnitude of the disorder 
potential at the TI surface using a self-consistent theory based on the Thomas-Fermi approximation 
for screening by the Dirac mode. Simple formulas are presented for the mean squared potential both 
at the Dirac point and far from it, as well as for the characteristic size of electron/hole puddles at 
the Dirac point and the total concentration of electrons/holes that they contain. We also derive an 
expression for the autocorrelation function for the potential at the surface and show that it has an 
unusually slow decay, which can be used to verify the bulk origin of disorder. The implications of 
our model for the electron conductivity of the surface are also presented. 



I. INTRODUCTION 

The three-dimensional (3D) topological insulator 
(TI^EHU has gapless surface states with a Dirac-like spec- 
trum, which host a number of interesting quantum trans- 
port phenomengP^ These surface states are influenced 
by the presence of a random Coulomb potential, which is 
believed to determine the mobility of surface electrons^. 
Recently, the random potential at the surface of typi- 
cal TIs was studied directly by spectroscopic mapping 
with a scanning tunneling microscope^. It was shown 
that near the Dirac energy random fluctuations of the 
potential have a Gaussian-like distribution with a width 
~ 20 - 40 meV. For a theoretical interpretation of their 
results, Ref. GO used a model of random charges with two- 
dimensional (2D) concentration situate d in a plane 
parallel to the surface at a distance d from it! 10 * 11 ! that is 
self-consistently screened by the electrons of the surface 
Dirac mode. Originally, this model was suggested to de- 
scribe disorder in graphene, where random charges can 
be assumed to be localized on the surface of a nearby sub- 
strate. It has since been extended to describe the surface 
of a 3D Tpl. 

In this paper we explore a different model of Coulomb 
disorder in TIs. We assume that the bulk of the TI 
is a completely compensated semiconductor with equal 
3D concentration N of donors and acceptors, which are 
randomly-distributed throughout the bulk . Thi s model 
is mathematically simpler than the 2D on e 1 10 * 11 ! because 
impurities are characterized by only one parameter, N, 
rather than two parameters, rij and d. An analysis of 
this 3D model is presented below, but before turning to 
it we would like to give arguments for the 3D model that 
are specific for known TIs. 

First, such a model is justified by current methods 
of preparation of TI crystals. Typically, as-grown TI 
crystals are heavily doped n-type semiconductors with 
N ~ 10 19 donors per cm 3 . The Fermi level of such a crys- 
tal is high in the conduction band. In order to bring the 
Fermi level down to the middle of the gap and increase 



the bulk resistivity, the TI is compensated by acceptors 
with concentration close to that of the donors, N. Below 
we assume that these donors and acceptors are randomly 
distributed in space. This is indeed usually the case for 
samples made by cooling from a melt, where the distri- 
bution of impurities in space is a snapshot of the dis- 
tribution at much higher temperature, when diffusion of 
impurities practically freezes^ 3 -'. In semiconductors with 
a narrow enough forbidden gap E gi at this temperature 
there is a concentration of intrinsic carriers larger than 
the concentration of impurities. These intrinsic carri- 
ers screen the Coulomb interaction between impurities, 
so that the impurities remain randomly distributed in 
space. As the melt is cooled to the point where intrin- 
sic carri ers re combine, the impurities are left in random 
positionPESl. If diffusion freezes at T ~ 1000 K it is 
reasonable to assume that impurities are randomly posi- 
tioned in a semiconductor with E g < 0.3 eV. 

Second, the model of a completely compensated TI 
with randomly positioned charges was recently tested by 
calculation of the activation energy A of its bulk resis- 
tivity and comparison to experiment)^. The standard 
expectation, which assumes flat valence and conduction 
bands, was that when the Fermi level is moved to the 
middle of the gap the bulk of the TI becomes a good in- 
sulator with activation energy A — E g /2. In reality, in 
the Bi2Se3 and Bi2Te3 families of TIs, where E g ~ 0.3 eV, 
the activation energy A was founcP^ to be frustratingly 
small, with A « 0.1hE g . This unexpectedly small activa- 
tion energy was shown to be explainable within the model 
of random 3D donor and acceptor charges. Specifically, 
it was shown by numerical simulation^ that A w Q.\bE g 
results from band bending by the potential created by 
random 3D Coulomb impurities. For these reasons we 
consider our model of 3D randomly situated donors and 
acceptors to be an appropriate description of the TI bulk. 

Our primary result for this model is an expression for 
the amplitude of fluctuations of the electric potential en- 
ergy, r, at the TI surface as a function of the chemical 
potential, fi, measured relative to the Dirac point. In 
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particular, for /i = we show below that 



r 2 = 



Q 4/3 



( M = 0). 



(1) 



Here, — e is the electron charge, k is the effective dielectric 
constant, and a = e 2 / nhv is the effective fine structure 
constant, where h is the reduced Planck constant and v is 
the Dirac velocity. This expression describes screening of 
the disorder potential via the formation of electron and 
hole puddles at the TI surface. The characteristic size of 
these puddles is given by 



TV' 1 / 3 
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(2) 



and the corresponding total number of electrons (or 
holes) per unit area in surface puddles is given by 



(3) 



Eqs. - (§ are derived below, along with results cor- 
responding to large /i. Our result for the mean squared 
potential, T 2 , is plotted in Fig.[l]as a function of /i. Below 
we also derive a simple relation for the autocorrelation 
function of the potential at the TI surface, which has an 
unusually slow decay and can be used to verify the bulk 
origin of disorder. 

In addition to describing the disorder potential, we cal- 
culate the corresponding electron conductivity a of the 
surface, and we show that when the average electron con- 
centration n satisfies n S> n p , the conductivity is given 
by 



2^ 



,3/2 



h Q 2 ln(l/a) N 



(4) 



where e 2 /h is the conductance quantum. At much 
smaller electron concentrations, n <C n p , the conductiv- 
ity saturates at a value <r m j n , which we estimate as 



e 2 1 

h 7ra ln(l/a) 



(5) 



The remainder of this paper is organized as follows. In 
Sec. [Tl] we develop our self-consistent theory to describe 
the screened disorder potential at the TI surface. In Sec. 
injthese analytical results are compared to results from a 
numerical simulation of the TI surface. Sec. |IV| presents 
the implications of our model for the electron conductiv- 
ity of the surface. We conclude in Sec. |V] by comparing 
our theory with the recent experiments of Ref. 9j and by 
discussing the major assumptions of our theory and its 
implications for future experiments. 



II. SELF-CONSISTENT THEORY OF THE 
SURFACE DISORDER POTENTIAL 

In the limit where the potential varies slowly compared 
to the characteristic Fermi wavelength of electrons at the 



surface, the electric potential </>(r) can be described using 
the Thomas-Fermi (TF) approach: 



H = E f [n{r) 



(r). 



(6) 



Here, Ef(n) = hv-\/4ir\n\ sgn(n) = 

(e 2 /an)y/ 47rjnj sgn(n) is the local Fermi energy 
and n(r) is the 2D electron concentration at the point 
r on the surface. The TF approximation is justified 
whenever a <C 1, as we show below. In TIs such small 
a can be seen as the result of the large bulk dielectric 
constant Kb > 30. We note here that for describing 
the properties of the surface state, which exists at a 
dielectric discontinuity, one should use for the effective 
dielectric constant k the arithmetic mean of the internal 
and external dielectric constants. If the TI is in vacuum, 
then k — («b + l)/2 ~ Kb/2. 

When the chemical potential is large enough in magni- 
tude that /i 2 e 2 ((j> 2 ), where (...) denotes averaging over 
the TI surface, the relation Ef(n) can be linearized to 
read Ef[n(r)} ~ /j,+Sn(r)/v(fi). Here Sn(r) = n(r)—n,Q is 
the difference in the electron concentration relative to the 
state with zero electric potential, n = a 2 k 2 /i 2 / (Aire 4 ) , 
and v(f-i) = a 2 K 2 \/j,\/(2ire 4 ) is the density of states at 
Ej = [i. From this density of states one can define a 
screening radius r s — K/2ire 2 v = e 2 /a k/j, that char- 
acterizes the distance over which fluctuations in the 
Coulomb potential are screened by the surface. The 
TF approximation is valid when the Fermi wavelength 
A 



f ~ n Q ~ e 2 /aK/i is much smaller than r s , which 
gives the condition a«l. 

One can understand qualitatively the magnitude of 
the potential fluctuations, T, using the following sim- 
ple argument. For a given point on the TI surface, one 
can say that only impurities within a distance R < r s 
contribute to the potential; those impurities at a dis- 
tance R 3> r s are effectively screened out (one can say 
that they are screened by their image charges in the 
"metallic" TI surface). Impurities with R < r s , on 
the other hand, are essentially unscreened. There are 
~ Nr^ such impurities, and their net charge is of or- 
der Q ~ e\J iVY 3 , with a random sign. The absolute 
value of the potential at the surface is then ~ Q / nr s , so 
that T ~ eQ/nr s ~ (e 2 iV 1 / 3 /K)(iVr 3 ) 1 / 6 - yJe 2 N/nv ~ 
Throughout this paper we assume that 
all donors and acceptors in the bulk are ionized, and we 
ignore the possible effects of bulk screening. This as- 
sumption is generally justified as long as the bulk chem- 
ical potential resides in the band gap, as we discuss in 
Sec. El 

In order to more accurately derive the value of T, one 
can start by considering the potential created by a single 
impurity charge +e. When such an impurity charge is 
placed a distance z from the TI surface (say, above the 
origin), it creates a potential <pi(r; z) that within the TF 
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approximation is given 



e exp -gz 
h(r;z) = - 



k J 1 

o 



7 TT J o(9 r ) d 9> 



(7) 



where Jo {%) is the zeroth order Bessel function of the first 
kind. At large z/r s , Eq. ^ can be expanded to give 



01 (r; z) 



k (r 2 + z 2 )3/2 ' 



(8) 



A simple physical derivation of Eq. ^ is based on the 
notion^ that for a distant impurity, such that z ^> r s , a 
surface with screening radius r s effectively plays the role 
of a metallic surface positioned below the real surface at a 
distance z — —r s /2. Equation (|Sj) can then be viewed as 
the sum of the potentials created by the original charge 
at a distance z above the plane and its opposite image 
charge at a distance z + r s below the plane, expanded to 
lowest order in r s /z. 

The total potential at the origin is (f)(0) = 
z i)i where the index i labels all impurity 
charges, g, is the sign of impurity i, and and Zi are the 
radial and azimuthal coordinates of its position. Under 
the assumption that all impurity positions are uncorre- 
cted and randomly-distributed throughout the bulk of 
the TI, the average of <fr 2 is given by 



l (r';z')] 2 2Nd 2 r'dz'. 



(9) 



Here, the quantity 2Nd 2 r'dz' describes the probability 
that the volume element d 2 r'dz' contains an impurity 
charge, and the integration is taken over the semi-infinite 
volume of the bulk of the TI. The width of the disorder 
potential at the TI surface, T, is defined by T 2 — e 2 ((f> 2 ). 
Inserting Eq. Q into Eq. ([9| and taking the integral then 



gives 



r 2 = 



e 2 N 2Tre i N 



a 2 K 3 \fi\ ' 



e 2 N l/3 



fi'O 



2/3 



(10) 



Eq. ([lOj) is correct so long as the fluctuations in the 
Coulomb potential energy are small compared to the 
chemical potential, or T <C this gives the condition 
written in parentheses. 

On the other hand, when is very small, the fluctua- 
tions in the Coulomb potential become large compared to 
the chemical potential, and one cannot talk about a con- 
stant density of states v or screening radius r s . Instead, 
the Fermi energy has strong spatial variations, and the 
random potential is screened by the formation of elec- 
tron and hole puddles at the surface. Nonetheless, one 
can define an average density of states (v) at the surface, 
which determines, self-consistently, the typical screening 
radius r s and the magnitude of the potential fluctuations 
at the TI surface. 

Consider, for example, the case [i = 0, where by sym- 
metry the average value of the potential (0) = 0. At 



any given point r on the surface, the potential 0(r) is 
the sum of contributions from many individual impurity 
charges, provided that the characteristic screening radius 
r s = K/2Tte 2 (v) > iV~ 1/3 . This implies that, by the 
central limit theorem, the value of the potential across 
the surface is Gaussian-distributed with some variance 
(0 2 ) = T 2 /e 2 that remains to be calculated. Within the 
TF approximation the local density of states at the point 
r is v\— e0(r)] = ea 2 K 2 |0(r)|/(27re 4 ), so that one can cal- 
culate the average density of states as 



v(-e(j>) 



exp [- 



J /2r 2 



a 2 K 2 Y 



v /2 7 rr 2 /e 2 



(P = o). 



(11) 



This result for (v) can be inserted into the first equality 
of Eq. (10), T 2 = e 2 N/n(v), to give a self-consistent re- 
lation for the amplitude of potential fluctuations^. This 
procedure gives the result first announced in the Intro- 
duction, Eq. ([I]). Substituting Eqs. ([I]) and (jTTj) into the 
expression for the screening radius, r s — K/2ire z (v) , gives 
Eq. §. 

One can also calculate the total concentration of elec- 
trons/holes in surface puddles, n p , implied by this result 
for r 2 . This is done by first inverting the TF relation, Eq. 

at n = to give n{4>) = (a 2 n 2 / Aire 2 )4> 2 sgn(0). Inte- 
grating this expression for n(<p) weighted by the Gaussian 
probability distribution for <f> gives 



n(<j>) 



exp [- 



2 /2r 2 



a 2 K 2 T 2 
87re 4 : 



v /27rT 2 7e 2 

(m = o). 



Substituting the result of Eq. |TJ for T 2 then gives Eq. 
(§. One can also combine this result for the residual 
electron/hole concentration, n p , with the expression for 
the screening radius, r s , to arrive at an estimate for the 
number of electrons/holes per puddle: M p ~ im p r 2 ~ 
7r/16a 2 . Apparently at small a puddles typically contain 
many electrons/holes, M p ^> 1. 

Our primarily results, outlined in Eqs. ([I]) - ([3]), are 
valid within the TF approximation so long as the typical 
Fermi wavelength, Xf ~ e 2 /anT, is much smaller than 
the typical screening radius, r s ~ e 2 /a 2 Kr, which again 
gives the condition a <C 1. Equations ([!]) and ^ were 
obtained in Ref. 21 without numerical coefficients within 
the context of graphene on a silicon oxide substrate. 

One can notice that our expressions for T 2 can be 
expressed more compactly by defining the dimension- 
less units r = T/Eq and ju = u/Eq, where E — 
e 2 N 1 / 3 /a 2 / 3 n. In these units Eqs. ([T]) and (10) can be 
written as 



T 2 = v^tt « 3.96, (p, = 0) 



(12) 
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and 



27T/I/ZI, (|M| » 1), 



(13) 



respectively, and the constant a does not enter explicitly. 
Eqs. ( 12 1 and ( 13 ) are plotted as the red dotted and 
dashed lines, respectively, in Fig. [I] One can similarly 
define a dimensionless screening radius r s = r s /r , where 
r = A" 1 / 3 /^ 3 , so that 



= 2~ 2/3 rs 0.63, (jl = 0) 



and 



I Ml 



m » i). 



(14) 



(15) 



At \i = 0, the screening radius r s describes the charac- 
teristic size of electron or hole puddles at the TI surface. 
More generally, r s plays the role of a length scale over 
which potential fluctuations at the surface are correlated. 
Such correlations can be discussed in a quantitative way 
by defining the potential auto-correlation function: 



C(r) = (0(R)0(R + r)) R , 



(16) 



where (...)r denotes averaging over the spatial coordi- 
nate R, and where by symmetry the correlation function 
depends on |r| = r only. Before proceeding to present 
numerical results, we first derive approximate analytical 
results for C(r), and show that spatial correlations in the 
potential have a n un usually slow decay. 

At r — 0, Eq. (16) reproduces the expression for (4> 2 ), 
so that C(0) = Tye 2 . At small enough distances that 
r <C r s , one can expect that the value of C(r) is deter- 
mined primarily by unscreened impurities that are within 
a distance r s from the surface, as explained above during 
the derivation of T 2 . On the other hand, at r > r s cor- 
relations are produced primarily by impurities that are 
relatively far from the surface, as can be seen from the 
following scaling argument. Consider two surface points 
separated by a distance r»r s . One can imagine draw- 
ing a cube of size r that extends into the bulk of the TI 
and which contains the two surface points on opposite 
edges of one of its faces. Such a cube contains ~ Nr 3 im- 
purities , and has a net impurity charge with magnitude 
q ~ e\/ Nr 3 and random sign. These impurity charges 
are located at a mean distance ~ r 3> r s above the sur- 
face and, therefore, b y Eq. Q , contribute a net potential 
~ qr s /nr 2 <~ {e/ n)yjNr 2 /r to both surface points. The 
square of this potential roughly gives the autocorrelation 
of the potential, C{r) ~ e 2 Nr 2 /n 2 r. 

A more careful expression for C(r) can be derived by 
writing 

C{r) = J 0i(r';2')<Mr' -r;z')2Nd 2 r'dz', (17) 

similar to Eq. Inserting the asymptotic expression 
of Eq. (pp for <f>i and evaluating the integral gives 

„, . 2ne 2 Nr 2 T 2 /e 2 

C(r )~ ^ = r/r.»l). 18 

K z r r/r s 



This result is plotted as the dashed line in Fig. [2] 

Eq. (181 implies an unusually slow decay of potential 



correlations at the surface, which, as explained above, 
arises from long-range fluctuations of the potential cre- 
ated by deep bulk impurities. This behavior can be con- 
trasted with the much faster decay of C(r) that would 
result from the 2D model of planar impurities^ C(r) ~ 
e 2 riidr 2 / n 2 r 3 . Thus, by studying C(r) experimentally 
by scanning tunneling microscopy, one can discriminate 
between disorder by bulk impurities and disorder by im- 
purities located in a layer close to the surface. 



III. NUMERIC SIMULATION 

So far we have presented analytical results for the mag- 
nitude of the potential at the surface and its autocor- 
relation function. These results are derived within the 
approximation of linear screening with a self-consistent 
screening radius r s . One can question the accuracy of this 
approach, particularly at small /x, where the density of 
states varies strongly from one point to another. There- 
fore, in order to verify the analytical results presented 
above, we implemented a simple simulation of a TI sur- 
face with adjacent point-like impurities and we solved 
numerically for the electric potential <f> at arbitrary Jl 
within the TF approximation. In these simulations, a 
square planar surface of dimension L x L is placed adja- 
cent to a volume of size L x L x L/2 with A^L 3 randomly- 
positioned impurities, each with a random sign. These 
impurities create a bare potential </> ex t(r) at the surface, 
and the self-consistent potential <f>(r) satisfies 



(r) =<f> ext (r)-[ d 2 v'^ 



, eSn(r') 



The TF equation, Eq. can be inverted to read 



6n(r) = 



Aire 3 



2^(r) + e[</)(r)] 2 sgn[0(r)], 



so that the self-consistent equation for the potential <f>(r) 
can be written 



<j>(r) = ^ext(r)- 



47re 2 



d 2 r' 



2/x0(r') + e[0(r')]' 



sgi#(r')]- 
(19) 



We find a solution </>(r) to Eq. (19 1 by dividing the 
surface into a discrete grid and using numerical iteration. 
Details of the iteration scheme, as well as our treatment 
of finite-size and finite-resolution effects, are given in the 
Appendix. Once we have obtained a numerical solution 
to 4>(r), the resulting variance of the disorder potential 
is calculated as 



r 2 = e 2 ((0-<0)) 2 ) 



(20) 



and the potential autocorrelation function is calculated 



using the definition in Eq. (16). All numerical results 
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presented below are calculated at a = 0.24 (as estimated 
for the experiments of Ref. |9]). Smaller a = 0.12 was also 
examined, and when presented in the dimensionless units 
of Eqs. |T2"| ) - ( [15] ) the results were identical to those of 
Figs. [I] and [2] to within our numerical error. 




2 4 6 

chemical potential, jl 

FIG. 1. (Color online) Variance in the disorder potential at 
the TI surface as a function of the chemical potential relative 
to the Dirac point. The dotted and dashed lines correspond to 
Eqs. ( 12 1 and ( 13 1, respectively. Open circles show the result 
of a numerical solution of Eq. ( |19[ ) ; error bars are smaller than 
the symbol size. 

In Fig. [I] the calculated value of V 2 is plotted as a 
function of Jl, along with the analytical asymptotic pre- 
dictions of Eqs. (fl2|) and (JT3J). (While Fig. [I] presents 



results only for positive Jl, results at Jl < are identi- 
cal due to electron-hole symmetry of the Dirac point.) 
The numerical results closely match the analytical the- 
ory at Jl = and at Jl 3> 1. One can noti ce, however, 
that at Jl s» Jl* = 2 2 / 3 , where Eqs. (12) and (13 1 become 
equal, r 2 develops a weak maximum. This maximum 
can be understood by considering that at Jl ~ Jl* the 
typical magnitude of the disorder potential, T, is similar 
to the typical Fermi energy, \x. As a result, screening is 
strongly asymmetric: positive fluctuations in potential, 
which increase the density of electrons, are screened more 
rapidly than negative fluctuations in potential, which de- 
plete the electron density and bring the system close to 
the Dirac point. The resulting distribution of the po- 
tential is skewed toward negative values of <f>, and this 
skewness produces a larger variance T 2 and a nonzero 
mean {(f)). As Jl is increased, of course, the width of the 
disorder potential becomes small relative to the typical 
Fermi energy, and screening becomes symmetric again. 

In Fig. [2] we plot the potential autocorrelation func- 
tion, C(r), as calculated from our simulation both at 
zero chemical potential, Jl — 0, and at large chemical po- 
tential, Jl = 8. In both cases the result compares well 



with Eq. (18) at large r/r s without the use of adjustable 
parameters. 
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Eq. (18) 
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FIG. 2. (Color online) The potential autocorrelation function, 
C(r), as a function of distance. Symbols correspond to data 
from our numeric simulation, both at fi = (squares) and 
at large /1 (triangles), while the dashed line is our analytical 
theory [Eq. |18l] for large r/r 3 . The ve rtical axis is scaled 
by the analytical result for T 2 [Eqs. ( 12 I and ( 13 1] and the 

[Eqs. 



horizontal axis is scaled by the analytical result for r s 
(14 1 and |15l] without fitting parameters. 



IV. CONDUCTIVITY 

In the previous sections we presented results for the 
disorder potential at the TI surface. In this section we 
discuss the implications of our 3D model for the conduc- 
tivity a of the surface. 

In the limit of large /x, where the electron density is 
only weakly modulated by the disorder potential, one can 
show using the Boltmzann kinetic equation that for elec- 
trons with a massless Dirac spectrum the conductivity is 
given b^l23H25J 



e 2 fir 



(21) 



Here e 2 /h is the conductance quantum and r is the mo- 
mentum relaxation time. In the limit of zero tempera- 
ture, the scattering rate 1 Jt can be found by integrating 
the squared scattering potential produced by a given im- 
purity over all impurities and over all scattering angles. 
More specifically, one can arrive at an expression for 1/t 
by taking the result for the scattering rate of a 2D layer 
of impurities with concentration rii at distance z [for ex- 
ample, Eq. (38) of Ref. replacing m with 2Ndz, and 
then integrating over all planes z containing impurities. 
This procedure gives 



1 kfOiK 

t Aithe 2 



2Ndz J dt) 




i 2 



<j)i(2k f sin-;z) 



(1 — cos" 



(22) 

In this equation, kf — an^/e 2 is the Fermi wave- 
length, 4>i{q]z) — (27re 2 /«g) exp[— qz]/{l + (qr s )~ 1 ] is the 
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screened potential (in momentum space) created by a sin- 
gle impurity at position z, and q = 2fc/sin(#/2) is the 
change in momentum associated with scattering by an 
angle 9. 



Evaluating the integral of Eq. ( 22 ) at small a gives 



1 e 2 N 
_^ aln( l /a) _ 



(23) 



Inserting this result for r into Eq. (21) and substituting 
kf 



'Aim yields the result for conductivity announced 
in the Introduction, Eq. Q. This expression can also be 
written in terms of the (dimensionless) chemical potential 
as 



h 4ira ln(l/a) ' 



(24) 



where Jl = /i/^N 1 / 3 /a 2 / 3 n), as defined in Sec. [nj 

Equation Q can be contrasted with the w idely-used 
result for the 2D model of charge impuritie d 10 1 12 1 23 1 2 -^, 
for which the conductivity is linearly proportional to the 
electron density: a/(e 2 /h) ~ (l/a 2 )(n/rii). This dif- 
ference can be understood conceptually by noting that, 
for large angle scattering, only those impurities at a dis- 
tance smaller than the Fermi wavelength, Ay ~ n -1 / 2 , 
contribute significantly to scattering. One can therefore 
define, roughly speaking, an effective 2D concentration 
of scattering impurities as N\f ~ N/n 1 / 2 . Inserting 
N/n 1 / 2 for rii gives a oc (l/a 2 )(n 3 ^ 2 /N), similar to Eq. 
Q. The remaining factor l/ln(l/a) in Eq. g is re- 
lated to low-angle scattering by distant impurities with 
z Xf. So far we are unaware of any transport data for 
TIs that shows a oc ?i 3 / 2 . Recent conductivity measure- 
ments on ultra-thin TIs (with thickness ~ 10 nm <C Xf) 
suggest 2 - a oc n, consistent with the 2D model of impu- 
rities. 

Our 3D model also produces a distinct result for the 
minimum conductivity cr m i n that appears in the limit of 
small average electron concentration. At small enough 
chemical potential that Jl <C 1, the surface breaks into 
electron and hole puddles, and one can think that the 
effective carrier concentration saturates at ~ n p [see Eq. 
([3])]. An estimate of cr m i n can therefore be obtained by 
setting ju = ju* = 2 2 / 3 in Eq. ( 24 ) , which gives the result 
of Eq. ([5]). 2D models of disorder impurities also pro- 
duce a minimum conductivity that is independent of the 
impurity concentration, but which has a different depen- 
dence on a. Specifically, at small a such models giv e ! 10 * 11 ! 
Cmin ~ (e 2 //i) ln(l/a). Our model suggests a minimum 
conductivity that is larger by a factor ~ [a In (l/a)] -1 . 



V. DISCUSSION 

In Sees, [n] and |III| we derived analytical expressions 
for the magnitude of the disorder potential, the screening 
radius, and the autocorrelation function, and we showed 



that these were consistent with numerical simulations. 
We now discuss the magnitude of T and r s implied by 
these expressions for typical TIs, which generally have an 
impurity concentration N ~ 10 19 cm' 3 . Typical values 
of the Dirac velocity and fine structure constant for TIs 
can be taken from Ref. 9, which reports hv = 1.3 eVA 
and estimates a = 0.24, which corresponds to k ~ 50 
(in agreement with infrared measurements on Bi2Se3, for 
example, which yielcP^ Kb ~ 100). Using these parame- 



W 3 / 



h'Ci 



2/3 



ters gives for our unit of energy Eq 
20meV, and r = N-^ 3 /a 4 / 3 ~ 30 nm. Thus, Eqs. ([12) 
and (14 1 imply T ~ 30meV and r s ~ 20 nm at the Dirac 
point, fi = 0. At large \fi\ > 30meV, both T 2 and r s 
decay as 

Throughout this paper we have worked within the as- 
sumption that bulk impurities are completely ionized, or 
in other words that there is no screening by conduction 
band electrons or valence band holes in the bulk. Such 
an assumption is valid when the chemical potential re- 
sides in the middle of a large bulk band gap. In this 
case donors or acceptors can only be neutralized by very 
large band bending!^, which brings the bulk conduction 
or valence band edges to the Fermi level (see, for exam- 
ple, Fig. 1 of Ref. ITFj) . Such fluctuations take place over 
a long length scale R g that scales as the square of the 
distance between the Fermi level and the nearest band 
edge. For example, if the Fermi level is in the center of 
the band gap, then R g = E 2 n 2 /8nNe 4 , which is on the 
order of hundreds of nanometers for typical TIs^. On 
the other hand, near the surface of the TI the potential 
fluctuations are screened much more effectively and over 
a much shorter distance, r s , by the (ungapped) surface 
states. As shown above, r s is typically < 20 nm, and the 
amplitude of surface potential fluctuations r ~ 30 meV 
<C Eg ~ 300 meV. One can therefore safely assume that 
near the surface there is no large band bending and one 
can indeed treat bulk impurities as completely ionized. 
The effect of bulk screening should appear only in the 
long-range behavior of the correlation function, r 3> R g , 
where the 1/r decay of C(r) is truncated and, as one can 
show, is replaced with C(r) ~ e 2 NR g r 2 / 'n 2 r 2 . 

We now compare our results for T and r s to the recent 
experiments of Ref. 9, Those authors examined sam- 
ples of doped Bi 2 Te3 and Bi 2 Se3 for which \x ~ 100 meV, 
and they found that the electric potential at the surface 
was well-characterized by a Gaussian distribution with a 
width ~ 20 - 40meV, which corresponds to T ~ 10 - 
20 meV. The characteristic length scale of potential fluc- 
tuations was measured as r s ~ 20 - 30 nm. Using the 
estimate above for Eq suggests that for these samples 
Jl ~ 6. Equations (13) and (15) then give T ~ 18meV 
and r, 



5 nm, which is in reasonable agreement with 
experiment. Further, evidence from Ref . |H] suggests that 
the surface disorder potential originates primarily from 
impurities deep below the TI surface, and that its magni- 
tude is relatively independent of the type of (monovalent) 
impurity present. These findings are again consistent 
with the theory presented here. So far, measurements 



7 



of the potential autocorrelation function C(r) have not 
been reported, but they can in principle be extracted 
from the measurements of Ref. [§] and compared to our 
prediction above. 

It is also worth mentioning that the theory we have 
presented here can be applied to graphene on a substrate 
having large dielectric constant (so that a = e 2 / nhv <C 1) 
and embedded bulk impurity charges. For this applica- 
tion one should only replace v everywhere with Av 1 since 
graphene 's spin and valley degeneracy give it a four times 
larger density of states at a given energy. This substitu- 
tion produces values of V 2 and r s for graphene that are 
four times smaller than what is written in Sees. Ill and fill 

Finally, we note that our theory ignores the possibility 
of screening by material outside the TI. For example, if 
the TI is placed next to a metal electrode or an ionic 
liquid^, then this external material can screen the large 
potential fluctuations created by the bulk, thereby de- 
creasing r and r s . 
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Appendix A: Numerical solution of the simulated 
self-consistent potential 



than setting </>W(r) = </>ncw(r) directly, we use a standard 
underrelaxation scheme with a damping parameter 7 to 
improve convergence of the solution: 



In Sec. |III| we presented results from a numeric solution 
of the potential at a simulated TI surface. In this Ap- 
pendix we discuss the details of our numerical method, 
including our iteration scheme and our treatment of finite 
size and finite resolution effects. 

In our simulation, NL 3 impurities, each with a ran- 
dom sign and random position, are placed in a volume 
with dimensions L x L x L/2 and open boundary condi- 
tions. One of the square faces of this volume is divided 
into a square grid with p(L + l) 2 grid points, where p is 
the grid resolution. In order to solve numerically for the 
potential at this surface for a given value of the parame- 
ters p, L, and p, we use a numerical iteration scheme that 
makes successive approximations ^"'(r) for the potential 
at each grid point r using Eq. (19 1. The first approxima- 
tion is made by evaluating the bare potential </> ext (r) at 
each grid point r, which we equate with (f>(°'(r). We then 
evaluate the right hand side of Eq. ( [l~9| for each r, which 
we denote 0oew(r), by inserting </>(°)(r) for </>(r). Rather 



^" +1 )(r)= 7 0(")(r) + (l- 7 )^(r). 



(Al) 



This process is continued iteratively, with </4ew(r) evalu- 
ated at each iteration by inserting <f, 
and used to create a revised estimate 



(™J(r) into Eq. ([19} 
(j)( n + 1 ) (r) according 
to Eq. (I Alh . The iteration is halted when the value of T 2 



associated with <p( n >(r) [see Eq. (20)] has converged to 
within 0.01%. For each value of the simulation parame- 
ters, all results are averaged over 100 random placements 
of the bulk impurity charges. Results presented above use 
7 = 0.5, but we verified that our convergent solution is 
independent of the value of 7 chosen for 0.2 < 7 < 0.98. 

Because of the long-ranged nature of the potential fluc- 
tuations created by bulk impurities, finite size effects in 
these simulations can be significant. The ideal numeri- 
cal result corresponds to the limit of an infinitely large 
simulation volume with an infinitely well-resolved spa- 
tial grid, L — > 00 and p — > 00. In practice, approaching 
this limit very closely can require an unrealistically large 
simulation. We therefore make use of an extrapolation 
method to estimate the value of T 2 corresponding to the 
L — > 00 and p — > 00 limit. Specifically, we find that for 
a given value of the chemical potential /i and grid reso- 
lution p, the variance in the potential can be well fitted 
to the equation T 2 (L) = T 2 ^ — A/L, where A is some 
positive constant. This dependence can be justified the- 
oretically by considering that the perimeter of the sim- 
ulated surface, which contains a fraction cx 1/L of the 
total number of grid points, has a smaller amplitude of 
the potential than the center of the grid by virtue of be- 
ing at the edge of the simulation volume. The value of 
r 2 ^ is then extracted by making a best fit to T 2 (L) us- 
ing a range of simulated sizes. Typically, our simulations 
use LN 1 / 3 = 20, 25, 30, 40, 50, 60 and we see a coefficient 
of determination R 2 > 0.95 for the linear fit of T 2 as a 
function of 1/L. (For a = 0.24 the theoretical screening 
radius at fi = is r s ps 4.2A~ 1 / 3 .) 

A similar extrapolation is also performed to evaluate 
the limit p — > 00. The (l/L)-extrapolated values of T 2 at 
a given grid density p are fitted to a linear function of 1/p, 
and the final estimate of T 2 for a given p is equated with 
the y-intercept of the corresponding line. The result of 
both of these extrapolations produces an estimate for T 2 
that is at most 11% different from the ones taken directly 
from our largest simulated sizes. Larger T 2 generally 
corresponds to larger extrapolation. 

The correlation function C(r) plotted in Fig. [2] is not 
the result of an extrapolation, but is calculated directly 
from a single set of simulations with large L sa 50r s and 
p m (4/r s ) 2 , averaged over 100 random placements of the 
bulk impurity charges for each value of p. 
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